Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/03b_fun_calib_ntd.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())Cost benefit analysis of the federal interventions that aim to reduce the harms of opioid crisis
Final Project | EPIB676
Github repo link: https://github.com/mariamelsheikh/cba_opioiduse
1 Background
Opioids are drugs used for pain relief and management but they can also be used to induce euphoria (ie feeling high). The epidemic of opioid-related overdose deaths and related harms began in Canada in 2016 and is still ongoing. This epidemic resulted in over 30,000 opioid-related overdose deaths since 2016, and in a huge burden on the healthcare system, just directly from opioid-related poisoning there has been over 34 000 opioid-related (excluding Quebec) and over 28 000 emergency medical services responses to suspected opioid-related overdose 2022 alone. It is important to also recognize that the harms of opioids epidemic are more than just opioid-related overdose deaths. (Federal 2023) This crisis is very complex with many factors affecting it but two of the main contributors to it are 1) illegal opioids: since the illicit drug market have become contaminated with strong opioids such as fentanyl (will be referred to them as illicit opioids) and 2) prescription opioids are also contributing to the crisis since they can be misused.(Health Canada 2020) As a result, through the Substance Use and Addictions Program (SUAP) and other related programs, the government is funding projects that support people who use drugs and aim to mitigate the harms and deaths resulting from opioid overdose epidemic based on scientific evidence Canada (2020).
Three of these interventions are:
- Increased access to overdose education and naloxone dispensation (Naloxone): Naloxone is a drug that reverses an opioid overdose and is proven to be safe. This intervention aims to distribute 58 000 naloxone kits across Canada and provide training and awareness on opioid-related overdose response (funding for it is $20 million dollars).(Canada 2021a)
- Opioid prescription guidelines: In 2017 stricter prescription guidelines were introduced for non-cancer chronic pain patients. These guidelines include 4 strong recommendations and 6 weak recommendations that aim to optimize non-opioid therapy and to restrict it to patients who do not have history of substance use.(Busse et al. 2017)
- Safer supply: which refers to services that aim to provide prescribed medications as a safer alternative to toxic illegal drug supply to people who are at high risk of overdose. This differs from agonist treatment (OAT) in terms of aims, population it serves and services. While OAT aims to treat substance use disorder and mitigates the effects of withdrawal by the use of different medications, Safer supply services aims to mitigate the harmful of effects of contaminated opioids from the illicit drug market, it is more flexible and does not aim to eventually stop drug use. (Canada 2021b)
This project aims to evaluate the future economic impact of these interventions by conducting a cost-benefit analysis to examine the difference in healthcare costs between the interventions and their projected effects on reducing opioid-related deaths and other important health outcomes. This project aims to provide evidence based on the currently available data on the effectiveness of these interventions.
2 Methods
This project is a team project, I worked on the modeling part but many of the decisions and the work upon which the modeling was built was done by other teams members (DP, MS, IW). The model was designed and conceptualized by DP. A systematic review was conducted by DP, MS and IW to extract the effects of the interventions from the literature. Transition probabilities, and healthcare costs were estimated by DP, either by extracting them from the literature or based on scientific assumptions in consultation with experts working in the field and/or people with lived experiences. I worked on the modeling part of the project with DP, and many of the decisions made for the modelling part were made in consultation with DP. Data sources for the parameters were either published journal papers, publicly available data and reports from different federal and provincial agencies, as well as other health organizations.
An open cohort Markov model that consists of 31 states representing the pathways of opioid use in Canada: “No risk of opioid-related overdose”, “At-risk of opioid-related overdose”, “Overdose and Sequelae” and “Post-brain injury opioid use states” was used for this analysis (fig. 1). The model runs over 15 years (2015 to 2029) with monthly cycles for population aged 15+ years (designed and conceptualized by DP). It is an open cohort model that takes into account the changes in population over time by adding people to the first state (Pain free, no use) every cycle (The yearly increase was assumed to be constant throughout the year, meaning that each increase was divided by 12. As for years beyond 2022, projections of population size in Canada were used (Government of Canada 2019) and proportion of those 15+ was estimated based on previous data, then a similar approach was used to estimate the monthly increase in the cohort).
Table 1: List of variables corresponding to the states in the model
Code
flextable(ini_states_tbl_org %>%
select(`Description `, `Variable `)) |>
autofit()Description | Variable |
|---|---|
Pain free, no use | BN_PN |
Acute pain, no use | BN_ACUTE |
Chronic (non-cancer) pain, no use | BN_CHRONIC |
Cancer, no use | BN_CANCER |
Other (i.e., cough, diarrhea, sickle cell disease), no use | BN_OTHER |
Palliative, no use | BN_PALLIATIVE |
Acute pain, Rx use | BPO_ACUTE |
Chronic (non-cancer) pain, Rx use | BPO_CHRONIC |
Cancer, Rx use | BPO_CANCER |
Other (i.e., cough, diarrhea, sickle cell pain), Rx use | BPO_OTHER |
Palliative, Rx use | BPO_PALLIATIVE |
Rx opioid misuse | BPO_MISUSE |
Illicit opioid use | BI_ILLICIT |
Detox / Withdrawal Management | BS_DETOX |
OAT initiation | BS_OAT_INI |
OAT maintenance | BS_OAT_MAINT |
OAT / Safe Supply | BS_OAT_SS |
Safe Supply | BS_SS |
Overdoes - illicit | BO_OD_ILLICIT |
Overdose - misuse | BO_OD_RX |
Moderate brain injury | BO_MOD_BI |
Severe brain injury | BO_SEVERE_BI |
Severe brain injury Out | BO_SEVERE_BI_OUT |
Opioid-related death | BO_OD_DEATH |
Death | BO_DEATH |
R: Illicit opioid use | BR_ILLICIT |
R: OAT initiation | BR_OAT_INI |
R: OAT maintenance | BR_OAT_MAINT |
R: OAT / Safe Supply | BR_OAT_SS |
R: Safe Supply | BR_SS |
R: Illicit opioid overdose | BR_OD_ILLICIT |
Characteristics of this model:
- One of the unique aspects of this model is that it distinguishes between chronic non-cancer pain, and cancer pain and it includes palliative care, which are important pathways that would capture the spillover effects of at least one of the interventions we are examining (prescription guidelines).
- We included brain injuries as part of the overdose sequelea. We did not distinguish between mild brain injury and no brain injury, both transition from overdose back into the pre-overdose/brain injury opioid-use states, however we distinguished between moderate and severe brain injury. For severe brain injury as a result of an opioid-related overdose we assumed that these individuals will not be part of the model anymore, they can either transition to death or transition to an “out” state where they can stay there until they die (we created the “out” state so that the high brain injury healthcare costs would not be incurred every cycle for those who did not die). For the moderate brain injury state, we created another set of repeat states for the opioid use states (repeat illicit opioid use, repeat OAT initiation, repeat OAT maintenance, repeat OAT/Safer supply, repeat Safer supply), since the transition probabilities would be different for those who experienced moderate brain injury.
- It distinguishes between opioid-related overdose deaths and deaths from any other causes (overall deaths).
- It distinguishes between prescription opioid misuse (defined as prescription-graded opioids used either by individuals other than those they were prescribed for, or used in a way other than prescribed) and illicit opioid use (opioids from illicit market), and distinguishes between overdose as a result of an prescription opioid and overdose from illicit opioids.
- It accounts for fentanyl contamination in the illicit drug market by including a multiplier of (1.005)^(t-36) (t-36 is number of month) to the probabilities of transitioning from illicit opioid use to illicit use overdose and transitioning from illicit use overdose to overdose death, from 2018 to 2020. Then beyond 2020 it remained constant.
2.1 Alternatives
This model was run for 5 scenarios: 1) No interventions, 2) Naloxone, 3) Prescription guidelines, 4) Safer Supply, and 5) all interventions combined. The effects of each of the interventions were estimated by other team members who conducted a systematic review of the literature (DP, MS, IW). Transition probabilities were changed for the naloxone and prescription guidelines interventions (and as a result the corresponding 1-rest transition probabilities will be changed).
For the naloxone intervention the following transition probabilities were changed starting from January 2017:
- illicit opioid overdose to overdose death was reduced by 5% (1-rest: illicit opioid overdose to illicit opioid use), and
- prescription opioid overdose to overdose deaths was reduced by 3% (1-rest: prescription opioid overdose to prescription opioid misuse).
For the prescription guidelines the following transition probabilities were changed starting from May 2017:
- pain free to acute pain, prescription use was reduced by 43% (1-rest: pain free stay in pain free),
- cancer no use to cancer prescription use was reduced by 55% (1-rest: cancer no use stay in cancer no use),
- acute pain prescription use to chronic pain, prescription use was reduced by 22% (1-rest: acute pain prescription use to pain free), and
- chronic pain prescription use to chronic pain no use was increased by 4% (1-rest: chronic pain prescription use stay in chronic pain prescription use).
For safer supply to account for the effect of the pilot programs introduced in 2016, 1000 people were moved from illicit drug use state to OAT/Safer supply state in January 2016. Then to account for compassionate prescribing due to the COVID-19 pandemic 10,000 people were moved from illicit drug use state to OAT/Safer supply state in April 2020. (Wyton 2022) Before January 2016 all transitioning into and out of safer supply states (Safer Supply, OAT/Safer Supply, Repeat OAT/Safer Supply and Repeat Safer Supply). Starting from January 2016, all transitioning out of OAT/Safer Supply (and staying in) were changed to their corresponding non-zero values. Until in 2023, when the intervention is supposed to be implemented, all transitioning probabilities into and out of safer supply states were changed to their non-zero values.
For the cost-benefit analysis the effects of these interventions individually and all together was of interest, therefore 4 different comparisons were made: 1) Naloxone vs No Interventions, 2) Safer Supply vs No Interventions, 3) Prescription Guidelines vs No Interventions and 4) All Interventions vs No Interventions.
2.2 Outcomes
Primary Outcomes:
- Cost-benefit analysis: Opioid-related costs from the healthcare perspectives for all states were used in this analysis. These include costs for hospitalization, ED visits, physician billing, paramedic services, drugs, and other healthcare system related costs. These were identified from publicly available data, as well as published systematic reviews, studies, and reports. All costs were discounted using 3% discount factor annually.
- Epidemiological measures: cumulative overall deaths and opioid-related overdose deaths.
Secondary Outcomes - more epidemiological measures:
- Weighted yearly prevalence of 1) severe brain injuries due to opioid overdose, 2) moderate brain injuries due to opioid overdose, 3) substance use treatment and 4) prescription opioid use.
For the primary outcomes changes and percent changes were calculated for the 4 different comparison groups while the secondary outcomes weighted yearly prevalence was calculated for each of th 5 scenarios.
2.3 Calibration and Uncertainty analysis
After designing the model structure, initial populations, and transitions, we calibrated the model using data from 2015 to 2021. Our calibration targets were total deaths (excluding opioid-related overdose deaths), the total number of people in OAT, and proportions of prescription opioid use. To obtain this data, we used various sources, including Statistics Canada, the Canadian Institute for Health Information’s report on Opioid prescribing in Canada, and scientific literature (DP).
Even though we had opioid-related overdose deaths numbers, we did not use them as targets because they are directly affected by different interventions that has been implemented so we could not calibrate the “no intervention” parameters to them. However, we used them to assess if the estimates of opioid-related overdose deaths from the model of the “no intervention” scenario are sound (they should be more than the targets).
We selected parameters to calibrate based on a one-way sensitivity analysis, where we changed one transition probability at a time, using lower and upper ranges (which were either estimated from the literature and if not available +/-25% range was use), and observed how this impacted predicted outcomes (overall deaths and overdose deaths). The one-way sensitivity analysis was conducted twice, 1) accounting for the fentanyl contamination (Scenario 1) and 2) without accounting for fentanyl contamination (Scenario 2). We included transition probabilities that resulted in a notable change in predicted outcomes (defined as 2.5%) in the set of parameters to be calibrated, in addition to all transition probabilities to overall deaths or opioid-related overdose deaths (a total of 38 parameters for scenario 1 and 39 paramters for scenario 2). We then conducted a Maximum A Posteriori (MAP) estimation using Nelder Mead, a sampling method that aims to find the global optimum twice once with the model for both scenarios. This allowed us to select the parameter set that best approximated the targets and provided us with information to inform changes in the prior information for some of the parameters. Overall, we updated the following baseline transition probabilities based on MAP parameter set:
- Pain free to prescription opioids misuse changed from 0.00032511 to 0.00220011
- Chronic pain, prescription use to prescription overdose from 0.00097898 to 0.00160398
- Detox to illicit opioid use from 0.1818112 to 0.22181118
- Detox to OAT initiation from 0.80000000 to 0.82000000
- Palliative care, prescription use to death from 0.37500000 to 0.38250000
- OAT maintenance stay in OAT maintenance from 0.70000000 to 0.71000000
Then we used the Sample Importance Resampling (SIR) Bayesian calibration method after updating the prior information based on the MAP results, and we ran the model using the sample sets (that were randomly sampled from the priors defined for the parameters: either estimated from the literature or using a range of +/-5%), estimated likelihood for each set, and weighted the sets by likelihood to approximate posteriors.
Probabilistic sensitivity analysis (PSA) was conducted on primary outcomes to account for the uncertainty in our results. The SIR generated posteriors for the parameters that were calibrated were used as priors and for other transition probabilities, a uniform distribution was used with ranges either from the literature and if not available, +/-5% were used as ranges. From these priors, probabilistic samples were randomly sampled and the analysis was simulated 10 000 times.
3 Results
3.1 Uncalibrated No Intervention model
Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
trace_tbl_inc <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count")
trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_od_deaths
trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_mod_bi
trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_sev_bi
trace_tbl_inc$state <- factor(trace_tbl_inc$state,
levels = v_state_names,
labels = v_state_names)
p4 <- ggplot(trace_tbl_inc, aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl_inc %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal)Code
plotly::ggplotly(p4)Fig 2a: Trace plots for the states in the model
Code
p5 <- ggplot(trace_tbl_inc %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = c(2015:2029)),
aes(x = year, y = count)) +
geom_line() +
ylab("Total opioid-related overdose deaths") + xlab("Year")
# +
# ggtitle("Total opioid-related overdose deaths over time")
plotly::ggplotly(p5)Fig 2b: Cumulative Opioid-realted overdose deaths over 15 year horizon
The results from the uncalibrated model (for “no intervention”) show that over time with total opioid-related overdose deaths would increase from 1837 in 2015 to 117 649 by the end of 2029 and overall deaths would increase from 23 089 in 2015 to 4 517 061. Furthermore, illicit opioid use increased from 290 063 to 2 485 660, while prescription opioid misuse increased initially then plateaued and slightly decreased. Additionally, there was a decrease in the number of individuals in the chronic pain with no opioid use and with prescription opioid use.
3.2 Comparison between SIR sample, uncalibrated sample, uncalibrated point estimates and calibration targets
The model approximated the calibration targets well: even though the distribution of prevalence of opioid prescription did not include the target values, the difference between the values was negligible. The distribution for overall deaths included the target values and the distribution for opioid-related deaths was larger than the target values (what we wanted the model to result in). However the distribution of OAT counts does not include the target values and that might be either because these values are not accurate enough to represent the values estimated by the model or that the calibration process did not work as well. Overall, the estimates from the model seem to be approximating values found in the literature.
3.3 Results for the primary outcomes (PSA)
Code
m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
round(ci_cost_diff_nalox[1]/1000000, 0),", ",
round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
round(ci_cost_diff_ss[1]/1000000, 0),", ",
round(ci_cost_diff_ss[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
round(ci_cost_diff_pg[1]/1000000, 0),", ",
round(ci_cost_diff_pg[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
round(ci_cost_diff_all[1]/1000000, 0),", ",
round(ci_cost_diff_all[2]/1000000, 0), "]"))
mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
round(ci_death_diff_nalox[1], 0),", ",
round(ci_death_diff_nalox[2], 0), "]"),
paste0(round(ci_death_diff_ss[3], 0), " [",
round(ci_death_diff_ss[1], 0),", ",
round(ci_death_diff_ss[2], 0), "]"),
paste0(round(ci_death_diff_pg[3], 0), " [",
round(ci_death_diff_pg[1], 0),", ",
round(ci_death_diff_pg[2], 0), "]"),
paste0(round(ci_death_diff_all[3], 0), " [",
round(ci_death_diff_all[1], 0),", ",
round(ci_death_diff_all[2], 0), "]"))
mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
round(ci_oddeath_diff_nalox[1], 0),", ",
round(ci_oddeath_diff_nalox[2], 0), "]"),
paste0(round(ci_oddeath_diff_ss[3], 0), " [",
round(ci_oddeath_diff_ss[1], 0),", ",
round(ci_oddeath_diff_ss[2], 0), "]"),
paste0(round(ci_oddeath_diff_pg[3], 0), " [",
round(ci_oddeath_diff_pg[1], 0),", ",
round(ci_oddeath_diff_pg[2], 0), "]"),
paste0(round(ci_oddeath_diff_all[3], 0), " [",
round(ci_oddeath_diff_all[1], 0),", ",
round(ci_oddeath_diff_all[2], 0), "]"))
effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
mean_death_diff, mean_oddeath_diff)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
mean_death_diff = "Deaths ",
mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present | Deaths | Opioid-related Overdose Deaths |
Naloxone | 209 [107, 362] | 625 [414, 926] | -10089 [-15009, -6697] |
Safer Supply | 4216 [3651, 4821] | -14180 [-17532, -11074] | -3221 [-4105, -2529] |
Prescription Guidelines | -25455 [-31480, -19866] | 14417 [8610, 20565] | -5458 [-11045, -1624] |
All Interventions | -21058 [-27059, -15403] | 960 [-5687, 7753] | -18544 [-28496, -11309] |
Code
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
round(ci_cost_diff_per_nalox[1], 2),", ",
round(ci_cost_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_cost_diff_per_ss[3], 2), " [",
round(ci_cost_diff_per_ss[1], 2),", ",
round(ci_cost_diff_per_ss[2], 2), "]%"),
paste0(round(ci_cost_diff_per_pg[3], 2), " [",
round(ci_cost_diff_per_pg[1], 2),", ",
round(ci_cost_diff_per_pg[2], 2), "]%"),
paste0(round(ci_cost_diff_per_all[3], 2), " [",
round(ci_cost_diff_per_all[1], 2),", ",
round(ci_cost_diff_per_all[2], 2), "]%"))
mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
round(ci_death_diff_per_nalox[1], 2),", ",
round(ci_death_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_death_diff_per_ss[3], 2), " [",
round(ci_death_diff_per_ss[1], 2),", ",
round(ci_death_diff_per_ss[2], 2), "]%"),
paste0(round(ci_death_diff_per_pg[3], 2), " [",
round(ci_death_diff_per_pg[1], 2),", ",
round(ci_death_diff_per_pg[2], 2), "]%"),
paste0(round(ci_death_diff_per_all[3], 2), " [",
round(ci_death_diff_per_all[1], 2),", ",
round(ci_death_diff_per_all[2], 2), "]%"))
mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
round(ci_oddeath_diff_per_nalox[1], 2),", ",
round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
round(ci_oddeath_diff_per_ss[1], 2),", ",
round(ci_oddeath_diff_per_ss[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
round(ci_oddeath_diff_per_pg[1], 2),", ",
round(ci_oddeath_diff_per_pg[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
round(ci_oddeath_diff_per_all[1], 2),", ",
round(ci_oddeath_diff_per_all[2], 2), "]%"))
effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
mean_death_diff_per, mean_oddeath_diff_per)
effects_tbl_per |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff_per = "Discounted Net Present Costs",
mean_death_diff_per = "Deaths ",
mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Percent Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths | Opioid-related Overdose Deaths |
Naloxone | 0.03 [0.01, 0.05]% | 0.01 [0.01, 0.02]% | -3.46 [-4.13, -3.13]% |
Safer Supply | 0.54 [0.46, 0.61]% | -0.31 [-0.38, -0.24]% | -1.11 [-1.97, -0.66]% |
Prescription Guidelines | -3.23 [-3.94, -2.55]% | 0.32 [0.19, 0.45]% | -1.81 [-2.76, -0.93]% |
All Interventions | -2.67 [-3.39, -1.98]% | 0.02 [-0.13, 0.17]% | -6.43 [-7.16, -5.58]% |
Code
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all,
ci_death_diff_per_nalox, ci_death_diff_per_ss,
ci_death_diff_per_pg, ci_death_diff_per_all,
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3)) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
# geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16, 17, 18)) +
# scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
theme_minimal())Fig. 4a: Percent change between the interventions and no intervention in terms of 3 primary outcomes: cost, overall deaths and opioid-related overdose deaths
Code
tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)) %>%
full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)), by = "Intervention") %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
geom_point(aes(color = Intervention), size = 1) +
geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention), height = 0) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_hline(yintercept = 0, linewidth = 0.25) +
geom_vline(xintercept = 0, linewidth = 0.25) +
xlab("Change in Costs Compared to No Intervention (%)") +
ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
theme_minimal())Fig. 4b: Percent change in opioid-related overdose deaths vs percent change in costs for all interventions compared to no interventions
The results show that 3 interventions individually and combined reduced the total number of opioid-related deaths compared to no interventions, with all interventions combined resulting in -6.413539 [-5.58, -7.16] % difference compared to no intervention scenario. All interventions combined and prescription guidelines resulted in reduction in costs as well, however Naloxone and Safer Supply increased costs slightly. Implementing all interventions combined reduced both costs and opioid-related overdose deaths. Regarding overall deaths though, except for safer supply, all other scenarios either increased overall deaths or result in inconclusive results.
3.4 Results for the secondary outcomes (point estimates)
Code
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)
prev_fun <- function(mod, i, v_states){
prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("prop_",
str_sub(v_yrs_prev, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev, 3, 4)))))
for (j in 1:yrs_prev) {
yr <- (v_yrs_prev[1] - 1) + j
if (length(v_states) > 1)
{
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
} else {
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
}
tot_pop_uncalib_tbl[, j] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
# final table for monthly prevalence over 15 years
prop_uncalib_tbl <- prop_uncalib %>%
pivot_longer(1:15, names_to = "grp", values_to = "value") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:15, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop)) %>%
mutate(intervention = mod_names[i])
# table with weighted yearly prevalence
prop_uncalib_wei_mean <- prop_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(value, tot_pop_val)) %>%
ungroup() %>%
mutate(intervention = mod_names[i])
return (list(prop_uncalib_tbl = prop_uncalib_tbl,
prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}
# Prevalence of severe brain injury
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>%
bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)Weighted Yearly Prevalence of Severe Brain Injury
Code
ggplotly(ggplot() +
geom_point(data = prev_sevbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
theme_minimal())Fig. 5a: Weighted yearly prevalence of severe brain injury over 15 years
Code
# Prevalence of moderate Brain injury
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>%
bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)Weighted Yearly Prevalence of Moderate Brain Injury
Code
ggplotly(ggplot() +
geom_point(data = prev_modbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5b: Weighted yearly prevalence of moderate brain injury over 15 years
Code
# Prevalence of substance use treatment
v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
"BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
"BR_OAT_SS", "BR_SS")
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>%
bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)Weighted Yearly Prevalence of substance use treatment
Code
ggplotly(ggplot() +
geom_point(data = prev_tx_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5c: Weighted yearly prevalence of substance use treatment over 15 years
Code
# Prevalence of prescription opioid use
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target) %>%
rename(value = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>%
bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value) %>%
mutate(intervention = "Target"))Weighted Yearly Prevalence of prescription opioid use
Code
ggplotly(ggplot() +
geom_point(data = prev_rx_opd_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5d: Weighted yearly prevalence of prescription opioid use over 15 years
Code
# Prevalence of illicit use
v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicit)$prop_uncalib_wei_mean
prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>%
bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)Weighted Yearly Prevalence of illicit opioid use
Code
ggplotly(ggplot() +
geom_point(data = prev_illicituse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5e: Weighted yearly prevalence of illicit opioid use over 15 years
Code
# Prevalence of misuse
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>%
bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)Weighted Yearly Prevalence of prescription opioid misuse
Code
ggplotly(ggplot() +
geom_point(data = prev_rxmisuse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5f: Weighted yearly prevalence of prescription opioid misuse over 15 years
Code
# Prevalence of overdose
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_od)$prop_uncalib_wei_mean
prev_od_yr_tbl <- noint_prev_od_yr_tbl %>%
bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)Weighted Yearly Prevalence of opioid-related overdose
Code
ggplotly(ggplot() +
geom_point(data = prev_od_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())Fig. 5g: Weighted yearly prevalence of opioid-related overdose over 15 years
The results show that none of the interventions affected prevalence of brain injury except for Safer Supply. Safer supply also resulted in an increase of prevalence in substance use treatment (which is expected since it is part of substance use treatment). Prescription Guidelines and All interventions combned resulted in a reduction in prevalence of prescription opioid use and prescription opioid misuse. No notable differences were observed for prevalence of opioid illicit drug use and opioid-related overdose.
4 Discussion and next steps
This analysis shows the promising effects of these interventions in reducing the direct harms of the opioid crisis by reducing opioid-related overdose deaths. From a healthcare perspective, when all interventions combined are implemented they would reduce costs while saving lives. However, the results show that there might be negative spillover effects regarding overall deaths which are not as promising. We hypothesize that these negative spillover effects are affecting marginalized sub-populations such as those suffering from chronic non-cancer pain, cancer pain and those in palliative care, we are interested in further examining the pathways that led to that increase in overall deaths.
One of the strengths of this analysis is that the model distinguished between chronic non-cancer and cancer pain and included palliative care which allows us to examine mechanisms and pathways that have not been previously examined, once the data is available. It also included brain injuries as one of the often neglected sequelae of overdose. The uncalibrated model estimated values that were similar to target values (a strength of the model), however, it seemed that the calibration did not make a significant difference (a limitation). The calibration might not have worked due to non-identifiability problem since we have many parameters that we were trying to calibrate and few calibration targets (Alarid-Escudero et al. 2018). However, that is due to the lack of data availability. This has been the main limitation of this project, there is a lack of data for the different parameters, many of which were based on scientific assumptions and knowledge of experts in the field. This model does not take into account heterogeneity in gender, age,or models the dynamic aspect of demand and supply of both naloxone and drugs in the illicit market or distinguished between different illicit drugs (heroin vs fentanyl vs other opioids). As a next step, we are interested in 1) repeating the analysis for only Ontario since they have more data available so we can better calibrate the model parameters by including more calibration targets, 2) conduct a value of information analysis to identify the parameters that we need to collect data on, 3) use British Columbia data to try to further understand the pathways that led to increase in overall deaths and repeat the analysis for this province similar to Ontario and 4) conduct the same analysis with societal costs in addition to healthcare costs.
This analysis will contribute to the evidence on the potential effectiveness of these federal interventions and hints at possible negative spillover effects that need to be examined further. It also reflects the need for collecting better data to be able to evaluate the effectiveness of these interventions as they are being implemented and expanded.
5 Appendix 1: Tables of the parameters used in the analysis
5.1 Initial Population counts in each state
Code
flextable(ini_states_tbl_org %>%
select(`Variable `, basevalue)) |>
autofit()Variable | basevalue |
|---|---|
BN_PN | 12,329,264 |
BN_ACUTE | 77,121 |
BN_CHRONIC | 6,000,000 |
BN_CANCER | 163,129 |
BN_OTHER | 6,003,590 |
BN_PALLIATIVE | 4,003 |
BPO_ACUTE | 216,300 |
BPO_CHRONIC | 3,530,940 |
BPO_CANCER | 123,666 |
BPO_OTHER | 594,870 |
BPO_PALLIATIVE | 14,158 |
BPO_MISUSE | 274,426 |
BI_ILLICIT | 290,063 |
BS_DETOX | 7,500 |
BS_OAT_INI | 2,461 |
BS_OAT_MAINT | 72,573 |
BS_OAT_SS | 0 |
BS_SS | 0 |
BO_OD_ILLICIT | 0 |
BO_OD_RX | 0 |
BO_MOD_BI | 0 |
BO_SEVERE_BI | 0 |
BO_SEVERE_BI_OUT | 0 |
BO_OD_DEATH | 0 |
BO_DEATH | 0 |
BR_ILLICIT | 0 |
BR_OAT_INI | 0 |
BR_OAT_MAINT | 0 |
BR_OAT_SS | 0 |
BR_SS | 0 |
BR_OD_ILLICIT | 0 |
5.2 Population increase every month
Code
flextable(pop_increase_tbl %>%
mutate(Year = as.character(Year)) %>%
rename("Increase per month" = increase_per_month)) |>
autofit()Year | Increase per month |
|---|---|
2015 | 18,637 |
2016 | 27,882 |
2017 | 32,458 |
2018 | 38,812 |
2019 | 41,247 |
2020 | 31,231 |
2021 | 19,700 |
2022 | 54,525 |
2023 | 48,524 |
2024 | 51,161 |
2025 | 50,779 |
2026 | 50,018 |
2027 | 49,005 |
2028 | 48,006 |
2029 | 47,094 |
5.3 Transition probabilities
Code
flextable(trans_prob_tbl_new %>%
filter(basevalue != 0) %>%
rename("From - state" = var_from,
"To - state" = var_to,
Basevalue = basevalue)) |>
autofit()From - state | To - state | Basevalue |
|---|---|---|
BN_PN | BN_PN | 999.00000000 |
BN_PN | BN_ACUTE | 0.00165257 |
BN_PN | BN_CANCER | 0.00050000 |
BN_PN | BN_OTHER | 0.00069143 |
BN_PN | BN_PALLIATIVE | 0.00030770 |
BN_PN | BPO_ACUTE | 0.00110106 |
BN_PN | BPO_MISUSE | 0.00220011 |
BN_PN | BI_ILLICIT | 0.00002940 |
BN_PN | BO_DEATH | 0.00032167 |
BN_ACUTE | BN_PN | 999.00000000 |
BN_ACUTE | BN_CHRONIC | 0.60000000 |
BN_ACUTE | BPO_MISUSE | 0.00004680 |
BN_ACUTE | BI_ILLICIT | 0.00001040 |
BN_ACUTE | BO_DEATH | 0.00053201 |
BN_CHRONIC | BN_PN | 0.00010000 |
BN_CHRONIC | BN_CHRONIC | 999.00000000 |
BN_CHRONIC | BN_CANCER | 0.00042924 |
BN_CHRONIC | BN_PALLIATIVE | 0.00030770 |
BN_CHRONIC | BPO_CHRONIC | 0.01000000 |
BN_CHRONIC | BPO_MISUSE | 0.00728000 |
BN_CHRONIC | BI_ILLICIT | 0.00010000 |
BN_CHRONIC | BO_DEATH | 0.00040000 |
BN_CANCER | BN_PN | 0.00010000 |
BN_CANCER | BN_CHRONIC | 0.00003130 |
BN_CANCER | BN_CANCER | 999.00000000 |
BN_CANCER | BN_PALLIATIVE | 0.00066950 |
BN_CANCER | BPO_CANCER | 0.55000000 |
BN_CANCER | BO_DEATH | 0.01390000 |
BN_OTHER | BN_PN | 0.00053731 |
BN_OTHER | BN_OTHER | 999.00000000 |
BN_OTHER | BPO_OTHER | 0.00349866 |
BN_OTHER | BO_DEATH | 0.00060864 |
BN_PALLIATIVE | BN_PALLIATIVE | 0.18500000 |
BN_PALLIATIVE | BPO_PALLIATIVE | 0.44000000 |
BN_PALLIATIVE | BO_DEATH | 999.00000000 |
BPO_ACUTE | BN_PN | 999.00000000 |
BPO_ACUTE | BN_CHRONIC | 0.60000000 |
BPO_ACUTE | BPO_CHRONIC | 0.06200000 |
BPO_ACUTE | BPO_MISUSE | 0.00600000 |
BPO_ACUTE | BO_OD_RX | 0.00743025 |
BPO_ACUTE | BO_DEATH | 0.00053201 |
BPO_CHRONIC | BN_CHRONIC | 0.01750000 |
BPO_CHRONIC | BPO_CHRONIC | 999.00000000 |
BPO_CHRONIC | BPO_MISUSE | 0.00470310 |
BPO_CHRONIC | BO_OD_RX | 0.00160398 |
BPO_CHRONIC | BO_DEATH | 0.00040000 |
BPO_CANCER | BN_CANCER | 0.00874161 |
BPO_CANCER | BPO_CHRONIC | 0.00177217 |
BPO_CANCER | BPO_CANCER | 999.00000000 |
BPO_CANCER | BPO_PALLIATIVE | 0.00100425 |
BPO_CANCER | BPO_MISUSE | 0.00470310 |
BPO_CANCER | BO_OD_RX | 0.00097898 |
BPO_CANCER | BO_DEATH | 0.01390000 |
BPO_OTHER | BN_OTHER | 0.02825010 |
BPO_OTHER | BPO_OTHER | 999.00000000 |
BPO_OTHER | BPO_MISUSE | 0.00470310 |
BPO_OTHER | BO_OD_RX | 0.00097898 |
BPO_OTHER | BO_DEATH | 0.00040000 |
BPO_PALLIATIVE | BPO_PALLIATIVE | 999.00000000 |
BPO_PALLIATIVE | BO_DEATH | 0.38250000 |
BPO_MISUSE | BN_PN | 0.01000000 |
BPO_MISUSE | BPO_MISUSE | 999.00000000 |
BPO_MISUSE | BI_ILLICIT | 0.00426953 |
BPO_MISUSE | BS_DETOX | 0.00100000 |
BPO_MISUSE | BS_OAT_INI | 0.00100000 |
BPO_MISUSE | BO_OD_RX | 0.00097898 |
BPO_MISUSE | BO_DEATH | 0.00064333 |
BI_ILLICIT | BN_PN | 0.00001000 |
BI_ILLICIT | BN_CANCER | 0.00064386 |
BI_ILLICIT | BI_ILLICIT | 999.00000000 |
BI_ILLICIT | BS_DETOX | 0.01000000 |
BI_ILLICIT | BS_OAT_INI | 0.00700000 |
BI_ILLICIT | BO_OD_ILLICIT | 0.01320000 |
BI_ILLICIT | BO_DEATH | 0.00080801 |
BS_DETOX | BN_PN | 999.00000000 |
BS_DETOX | BI_ILLICIT | 0.22181118 |
BS_DETOX | BS_OAT_INI | 0.82000000 |
BS_DETOX | BO_OD_ILLICIT | 0.01431324 |
BS_DETOX | BO_DEATH | 0.00309346 |
BS_OAT_INI | BI_ILLICIT | 999.00000000 |
BS_OAT_INI | BS_OAT_MAINT | 0.91700000 |
BS_OAT_INI | BO_OD_ILLICIT | 0.04050000 |
BS_OAT_INI | BO_DEATH | 0.00150000 |
BS_OAT_MAINT | BN_PN | 0.00460000 |
BS_OAT_MAINT | BI_ILLICIT | 999.00000000 |
BS_OAT_MAINT | BS_OAT_MAINT | 0.71000000 |
BS_OAT_MAINT | BO_OD_ILLICIT | 0.00420000 |
BS_OAT_MAINT | BO_DEATH | 0.00038326 |
BS_OAT_SS | BS_OAT_SS | 999.00000000 |
BS_SS | BS_SS | 999.00000000 |
BO_OD_ILLICIT | BI_ILLICIT | 999.00000000 |
BO_OD_ILLICIT | BS_DETOX | 0.26800000 |
BO_OD_ILLICIT | BS_OAT_INI | 0.05871300 |
BO_OD_ILLICIT | BO_MOD_BI | 0.03000000 |
BO_OD_ILLICIT | BO_SEVERE_BI | 0.00086172 |
BO_OD_ILLICIT | BO_OD_DEATH | 0.01900000 |
BO_OD_RX | BPO_MISUSE | 999.00000000 |
BO_OD_RX | BS_DETOX | 0.26800000 |
BO_OD_RX | BS_OAT_INI | 0.05871300 |
BO_OD_RX | BO_MOD_BI | 0.03000000 |
BO_OD_RX | BO_SEVERE_BI | 0.00086172 |
BO_OD_RX | BO_OD_DEATH | 0.01380000 |
BO_MOD_BI | BO_DEATH | 0.00874161 |
BO_MOD_BI | BR_ILLICIT | 999.00000000 |
BO_MOD_BI | BR_OAT_INI | 0.15100000 |
BO_SEVERE_BI | BO_SEVERE_BI_OUT | 999.00000000 |
BO_SEVERE_BI | BO_DEATH | 0.67000000 |
BO_SEVERE_BI_OUT | BO_SEVERE_BI_OUT | 999.00000000 |
BO_SEVERE_BI_OUT | BO_DEATH | 0.00090000 |
BO_OD_DEATH | BO_OD_DEATH | 1.00000000 |
BO_DEATH | BO_DEATH | 1.00000000 |
BR_ILLICIT | BO_DEATH | 0.01740000 |
BR_ILLICIT | BR_ILLICIT | 999.00000000 |
BR_ILLICIT | BR_OAT_INI | 0.02712532 |
BR_ILLICIT | BR_OD_ILLICIT | 0.03650000 |
BR_OAT_INI | BO_DEATH | 0.01740000 |
BR_OAT_INI | BR_ILLICIT | 999.00000000 |
BR_OAT_INI | BR_OAT_MAINT | 0.84075000 |
BR_OAT_INI | BR_OD_ILLICIT | 0.01285868 |
BR_OAT_MAINT | BO_DEATH | 0.00109940 |
BR_OAT_MAINT | BR_ILLICIT | 999.00000000 |
BR_OAT_MAINT | BR_OAT_MAINT | 0.84075000 |
BR_OAT_MAINT | BR_OD_ILLICIT | 0.01820000 |
BR_OAT_SS | BR_OAT_SS | 999.00000000 |
BR_SS | BR_SS | 999.00000000 |
BR_OD_ILLICIT | BO_MOD_BI | 999.00000000 |
BR_OD_ILLICIT | BO_SEVERE_BI | 0.03000000 |
BR_OD_ILLICIT | BO_OD_DEATH | 0.07700000 |
5.4 Healthcare costs for all states used in the analysis
Code
flextable(costs_tbl) |>
autofit()Variable | Cost |
|---|---|
BN_PN | 0.00 |
BN_ACUTE | 200.34 |
BN_CHRONIC | 70.93 |
BN_CANCER | 580.18 |
BN_OTHER | 70.93 |
BN_PALLIATIVE | 5,405.14 |
BPO_ACUTE | 15,002.06 |
BPO_CHRONIC | 103.05 |
BPO_CANCER | 2,186.73 |
BPO_OTHER | 103.05 |
BPO_PALLIATIVE | 2,904.79 |
BPO_MISUSE | 912.27 |
BI_ILLICIT | 97.10 |
BS_DETOX | 4,550.00 |
BS_OAT_INI | 1,500.88 |
BS_OAT_MAINT | 1,393.88 |
BS_OAT_SS | 1,432.77 |
BS_SS | 1,780.69 |
BO_OD_RX | 974.48 |
BO_OD_ILLICIT | 974.48 |
BO_MOD_BI | 6,825.46 |
BO_SEVERE_BI | 46,749.90 |
BO_SEVERE_BI_OUT | 1,295.96 |
BO_OD_DEATH | 0.00 |
BO_DEATH | 0.00 |
BR_ILLICIT | 97.10 |
BR_OAT_INI | 1,500.88 |
BR_OAT_MAINT | 1,393.88 |
BR_OAT_SS | 1,432.77 |
BR_SS | 1,780.69 |
BR_OD_ILLICIT | 97.10 |
6 Appendix 2: One way sensitivity analysis for the uncalibrated no intervention model
Code
library(here)
library(plotly)
theme_set(theme_minimal())
source(here("02_scripts/01_fun_data.R"))
t_owsa_ini_pop <- read.csv(file = here("01_data/ows_tbl_ini_pop.csv")) %>% select(-X)
ows_tbl_prob <- read.csv(, file = here("01_data/ows_tbl_prob.csv")) %>% select(-X)7 Initial population
7.0.0.1 Costs
Code
t_owsa_ini_pop_costs <- t_owsa_ini_pop %>%
select(costs_low, costs_high, name)
t_owsa_ini_pop_costs$name <- reorder(t_owsa_ini_pop_costs$name, t_owsa_ini_pop$costs_range)
t_owsa_ini_pop_costs <- t_owsa_ini_pop_costs %>%
pivot_longer(1:2, names_to = "group", values_to = "costs") %>%
mutate(costs_new = (costs - mod_basecase$total_net_present_cost)) %>%
filter(costs_new != 0) %>%
mutate(group = ifelse(group == "costs_high", "+25%",
ifelse(group == "costs_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_costs %>% filter(abs(costs_new) >= 500000000),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))7.0.0.2 Deaths
Code
t_owsa_ini_pop_deaths <- t_owsa_ini_pop %>%
select(deaths_low, deaths_high, name)
t_owsa_ini_pop_deaths$name <- reorder(t_owsa_ini_pop_deaths$name, t_owsa_ini_pop$deaths_range)
t_owsa_ini_pop_deaths <- t_owsa_ini_pop_deaths %>%
pivot_longer(1:2, names_to = "group", values_to = "deaths") %>%
mutate(deaths_new = (deaths - mod_basecase$m_M[181, "BO_DEATH"])) %>%
filter(deaths_new != 0) %>%
mutate(group = ifelse(group == "deaths_high", "+25%",
ifelse(group == "deaths_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_deaths %>% filter(abs(deaths_new) >= 5000),
aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Deaths from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use",
"ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))7.0.0.3 OD Deaths
Code
t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop %>%
select(od_deaths_low, od_deaths_high, name)
t_owsa_ini_pop_oddeaths$name <- reorder(t_owsa_ini_pop_oddeaths$name, t_owsa_ini_pop$od_deaths_range)
t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop_oddeaths %>%
pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>%
mutate(oddeaths_new = (oddeaths - (mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths))) %>%
filter(oddeaths_new != 0) %>%
mutate(group = ifelse(group == "od_deaths_high", "+25%",
ifelse(group == "od_deaths_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_oddeaths %>% filter(abs(oddeaths_new) >= 100),
aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in opioid-related overdose deaths from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use",
"ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use",
"ini_pop_BS_OAT_MAINT" = "OAT maintenance")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))8 Transition probabilities
8.0.0.1 Costs
Code
ows_tbl_prob_costs <- ows_tbl_prob %>%
select(costs_low, costs_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_costs$name <- reorder(ows_tbl_prob_costs$name, ows_tbl_prob$cost_range)
ows_tbl_prob_costs <- ows_tbl_prob_costs %>%
pivot_longer(1:2, names_to = "group", values_to = "costs") %>%
mutate(costs_new = (costs - as.numeric(mod_basecase$total_net_present_cost)),
costs_new_per = (costs_new/as.numeric(mod_basecase$total_net_present_cost)) * 100) %>%
filter(costs_new != 0) %>%
mutate(group = ifelse(group == "costs_high", "+25%",
ifelse(group == "costs_low", "-25%", group)))
name_cost_1 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 0.5) %>% select(name)
name_cost_2 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 1.5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_1$name),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))Code
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_2$name),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))8.0.0.2 Deaths
Code
ows_tbl_prob_deaths <- ows_tbl_prob %>%
select(deaths_low, deaths_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_deaths$name <- reorder(ows_tbl_prob_deaths$name, ows_tbl_prob$deaths_range)
ows_tbl_prob_deaths <- ows_tbl_prob_deaths %>%
pivot_longer(1:2, names_to = "group", values_to = "deaths") %>%
mutate(deaths_new = (deaths - as.numeric(mod_basecase$m_M[181, "BO_DEATH"])),
deaths_new_per = (deaths_new/as.numeric(mod_basecase$m_M[181, "BO_DEATH"])) * 100) %>%
filter(deaths_new != 0) %>%
mutate(group = ifelse(group == "deaths_high", "+25%",
ifelse(group == "deaths_low", "-25%", group)))
name_death_1 <- ows_tbl_prob_deaths %>% filter(abs(deaths_new_per) >= 0.5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_deaths %>% filter(name %in% name_death_1$name),
aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in deaths from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))8.0.0.3 OD-deaths
Code
oddeaths_bc <- as.numeric(mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths)
ows_tbl_prob_oddeaths <- ows_tbl_prob %>%
select(od_deaths_low, od_deaths_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_oddeaths$name <- reorder(ows_tbl_prob_oddeaths$name, ows_tbl_prob$od_deaths_range)
ows_tbl_prob_oddeaths <- ows_tbl_prob_oddeaths %>%
pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>%
mutate(oddeaths_new = (oddeaths - oddeaths_bc),
oddeaths_new_per = (oddeaths_new/oddeaths_bc) * 100) %>%
filter(oddeaths_new != 0) %>%
mutate(group = ifelse(group == "od_deaths_high", "+25%",
ifelse(group == "od_deaths_low", "-25%", group)))
name_oddeath_1 <- ows_tbl_prob_oddeaths %>% filter(abs(oddeaths_new_per) >= 5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_oddeaths %>% filter(name %in% name_oddeath_1$name),
aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in opioid-related overdose death from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_ILLICIT" = "Illicit opioid overdose to death",
"p_BO_OD_ILLICIT_BI_ILLICIT" = "Illicit opioid use to Illicit opioid overdose",
"p_BO_OD_RX_BPO_MISUSE" = "Rx opioid misuse to Rx opioid overdose",
"p_BR_OD_ILLICIT_BR_ILLICIT" = "R: illicit opioid use to R: illicit opioid overdose",
"p_BO_OD_DEATH_BR_OD_ILLICIT" = "R: Illicit opioid overdose to death",
"p_BO_OD_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to Illicit opioid overdose",
"p_BI_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to illicit opioid use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))9 Appendix 3: Comparison between Uncalibrated model and Calibration Targets
Code
library(here)
library(RColorBrewer)
source(here("02_scripts/01_fun_data.R"))
calib_target_tbl <- read_excel(here("01_data/markov_model_parameters_preprior.xlsx"),
sheet = "calibration_targets")
theme_set(theme_minimal())9.1 Prevalence of Rx opioids use
Code
# prevalence of people on prescribed opioids
v_yrs_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrs_prev_rx)
prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("prop_opioids_rx_",
str_sub(v_yrs_prev_rx, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev_rx, 3, 4)))))
for (i in 1:yrs_prev_rx) {
yr <- v_yrs_prev_rx[1] + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev_rx,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "Target"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_tbl,
aes(x = year_mon, y = target_val,
color = factor(year), fill = factor(year))) +
geom_point(data = prop_opioids_rx_target_tbl,
aes(x = year_mon, y = target_val),
color = "red") +
geom_point(data = prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon),
aes(x = year_mon, y = target_val),
color = "blue")+
xlab("Year_Month") + ylab("Prevalence of prescription opioid use") +
scale_fill_discrete(name = "year")+
scale_color_discrete(name = "year")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)The red points represent the observed target proportions of Rx opioids from CIHI report, the blue points represent an annual average (weighted by population during that month), the other coloured points represent monthly prevalence predicted from the model
Code
p2 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_wei_mean,
aes(x = year, y = target_val, color = factor(grp),
fill = factor(grp))) +
xlab("Year") + ylab("Prevalence of prescription opioid use")+
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p2)9.2 Deaths and OD-Deaths
Code
################ Deaths
# number of deaths
v_yr1_deaths <- c(2016:2020)
yrs_deaths <- length(v_yr1_deaths)
num_deaths_uncalib <- rep(NA, yrs_deaths)
for (i in 1:yrs_deaths){
yr <- (v_yr1_deaths[1] - 1) + i
num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
v_yr1_oddeaths <- c(2016:2021)
yrs_oddeaths <- length(v_yr1_oddeaths)
num_od_deaths_uncalib <- rep(NA, yrs_oddeaths)
for (i in 1:yrs_oddeaths){
yr <- (v_yr1_oddeaths[1] - 1) + i
num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total opioid-related overdose deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total opioid-related overdose deaths")) %>%
mutate(year = c(v_yr1_deaths, v_yr1_oddeaths),
group = "Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
facet_wrap(~target, scales = "free") +
theme(legend.position="bottom")
plotly::ggplotly(p3)Code
plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = ""))Code
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total opioid-related overdose deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total opioid-related overdose deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = ""))9.3 OAT
Code
# number of oat
v_yr1_oat <- c(2018:2021)
yrs_oat <- length(v_yr1_oat)
num_oat_uncalib <- rep(NA, yrs_oat)
for (i in 1:yrs_oat){
yr <- (v_yr1_oat - 1) + i
num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% "total_oat") %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = v_yr1_oat,
group = "Model"))
p3 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
theme(legend.position="bottom")
plotly::ggplotly(p3)9.4 Trace plot
Code
# build-in color palette
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
trace_tbl_inc <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count")
trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_od_deaths
trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_mod_bi
trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_sev_bi
trace_tbl_inc$state <- factor(trace_tbl_inc$state,
levels = v_state_names,
labels = v_state_names)
p4 <- ggplot(trace_tbl_inc, aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl_inc %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal)
plotly::ggplotly(p4)Code
p5 <- ggplot(trace_tbl_inc %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = c(2015:2029)),
aes(x = year, y = count)) +
geom_line() +
ylab("Total opioid-related overdose deaths") + xlab("Year")
plotly::ggplotly(p5)10 Appendix 4: Comparison between MAP calibrated model, uncalibrated model and calibration Targets
10.1 Scenario 1
It takes into account fentanyl contamination
Code
library(here)
source(here("02_scripts/03a_fun_calib_td.R"))
theme_set(theme_minimal())
opt_params <- readRDS(here("01_data/map_calib_td.RDS"))Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets
Code
mod_opt_map <- mod_calib(as.numeric(opt_params$par))10.2 Prevalence of Rx opioids use
Code
# prevalence of people on prescribed opioids
v_yrlast_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrlast_prev_rx)
prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("prop_opioids_rx_",
str_sub(v_yrlast_prev_rx, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrlast_prev_rx, 3, 4)))))
for (i in 1:yrs_prev_rx) {
yr <- (v_yrlast_prev_rx[1] - 1) + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
prop_opioids_rx_calib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_calib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_calib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[6], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Uncalibrated Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "target")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Calibrated Model"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>%
mutate(group = "Uncalibrated Model") %>%
bind_rows(.,prop_opioids_rx_calib_tbl %>%
mutate(group = "Calibrated Model"))
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
mutate(grp = "Target") %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Uncalibrated Model")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Calibrated Model"))
prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
labels = c("Calibrated Model", "Uncalibrated Model"),
levels = c("Calibrated Model", "Uncalibrated Model"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
color = group, fill = group)) +
geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)10.2.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use
Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p2 <- ggplot() +
geom_point(data = wprop_opioids_rx_tbl,
aes(x = year_mon, y = target_val,
color = grp, fill = grp), size = 1.5) +
xlab("Year") + ylab("Prevalence of Rx opioid use")+
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
plotly::ggplotly(p2)10.3 Deaths and OD-Deaths
Code
################ Deaths
# number of deaths
num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
yr <- 2015 + i
num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
num_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
yr <- 2015 + i
num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
num_od_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total OD-deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_deaths_calib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Calibrated Model"))
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p3)Code
p4 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p4)10.4 OAT
Code
# number of oat
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
yr <- 2017 + i
num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
num_oat_calib[i] <- sum(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_oat")) %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_oat_calib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Calibrated Model"))
oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p5 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total OAT counts") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p5)Code
cbind.data.frame(parameter = calib_param_tbl$var_params,
basevalue = calib_param_tbl$basevalue,
opt_map_nm = opt_params$par) %>%
filter(basevalue != opt_map_nm) parameter basevalue opt_map_nm
1 p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2 p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3 p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_DETOX__BI_ILLICIT 0.18181118 0.22181118
5 p__BS_DETOX__BS_OAT_INI 0.80000000 0.82000000
6 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.7100000010.5 Scenario 2
It does not take into account fentanyl contamination
Code
library(here)
source(here("02_scripts/03b_fun_calib_ntd.R"))
theme_set(theme_minimal())
opt_params_notimedepend <- readRDS(here("01_data/map_calib_ntd.RDS"))Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets
Code
mod_opt_map_ntd <- mod_calib_ntd(as.numeric(opt_params_notimedepend$par))10.6 Prevalence of Rx opioids use
Code
# prevalence of people on prescribed opioids
prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("prop_opioids_rx_",
str_sub(2015:2018, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("tot_pop_",
str_sub(2015:2018, 3, 4)))))
for (i in 1:length(2015:2018)) {
yr <- 2014 + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
prop_opioids_rx_calib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_calib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_calib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[6], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Uncalibrated Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "target")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Calibrated Model"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>%
mutate(group = "Uncalibrated Model") %>%
bind_rows(.,prop_opioids_rx_calib_tbl %>%
mutate(group = "Calibrated Model"))
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
mutate(grp = "Target") %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Uncalibrated Model")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Calibrated Model"))
prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
labels = c("Calibrated Model", "Uncalibrated Model"),
levels = c("Calibrated Model", "Uncalibrated Model"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
color = group, fill = group)) +
geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 2) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)10.6.0.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use
Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p2 <- ggplot() +
geom_point(data = wprop_opioids_rx_tbl,
aes(x = year_mon, y = target_val,
color = grp, fill = factor(grp)), size = 2) +
xlab("Year") + ylab("Prevalence of Rx opioid use") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p2)10.7 Deaths and OD-Deaths
Code
################ Deaths
# number of deaths
num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
yr <- 2015 + i
num_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
num_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
yr <- 2015 + i
num_od_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
num_od_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total OD-deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_deaths_calib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Calibrated Model"))
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p3)Code
p4 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p4)10.8 OAT
Code
# number of oat
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
yr <- 2017 + i
num_oat_uncalib[i] <- sum(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
num_oat_calib[i] <- sum(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_oat")) %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_oat_calib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Calibrated Model"))
oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p5 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total OAT counts") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p5)Code
cbind.data.frame(parameter = calib_param_tbl_ntd$var_params,
basevalue = calib_param_tbl_ntd$basevalue,
opt_map = opt_params_notimedepend$par) %>%
filter(basevalue != opt_map) parameter basevalue opt_map
1 p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2 p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3 p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.7700000011 Appendix 5: Comparison between SIR calibrated model, uncalibrated model and calibration targets
Code
library(here)
library(plotly)
source(here("02_scripts/01_fun_data.R"))
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("03_outputs/03_map_calibmod_vs_targets_td.R"))
calib_samples <- readRDS(here("01_data/calib_samples_sir.RDS"))
uncalib_samples <- readRDS(here("01_data/uncalib_sample_sir.RDS"))
opt_params_sir_prev <- readRDS(here("01_data/prev_outcome_data_sir.RDS"))
opt_params_sir_deaths <- readRDS(here("01_data/death_outcome_data_sir.RDS"))
opt_params_sir_oddeaths <- readRDS(here("01_data/oddeath_outcome_data_sir.RDS"))
opt_params_sir_oat <- readRDS(here("01_data/oat_outcome_data_sir.RDS"))
params_sir_prev_uncalib <- readRDS(here("01_data/prev_outcome_data_sir_uncalib.RDS"))
params_sir_deaths_uncalib <- readRDS(here("01_data/death_outcome_data_sir_uncalib.RDS"))
params_sir_oddeaths_uncalib <- readRDS(here("01_data/oddeath_outcome_data_sir_uncalib.RDS"))
params_sir_oat_uncalib <- readRDS(here("01_data/oat_sir_uncalib.RDS"))SIR calibration sample: first i sampled from prior 10000 sample sets, then i used them to calculated outcomes and likelihoods for each sample set, then i resampled from those 10000 with replacement with likelihood weights 10000 samples (unique sets = 6398), and ran the model on all of them and generated distributions for prevalence of rx opioid use, deaths, od deaths, total oat —- likelihood included deaths (gamma distribution), prevalence of rx opioid use (normal distribution instead of binomial), and total OAT
Uncalibrated sample: Here i resampled from the 10 000 sampled from the prior with replacement with equal weights 10000 sampled (unique sets = 6296), then ran the model on all of them and generated distributions for prevalence of rx opioid use, deaths, od deaths and total oat
11.1 Prevalence of Rx opioid use
Code
dt_prev <- as.data.frame(params_sir_prev_uncalib) %>%
pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_prev) %>%
pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
prev_year = c("yr15", "yr16", "yr17", "yr18")
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
rename(value = target_val) %>%
mutate(grp = "Target",
prev_year = ifelse(grepl("15", year_mon),2015,
ifelse(grepl("16", year_mon), 2016,
ifelse(grepl("17", year_mon), 2017,
ifelse(grepl("18", year_mon), 2018, year_mon)))),
prev_year = factor(prev_year)) %>%
select(-year_mon) %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(prev_year = factor(year),
grp = "Uncalibrated Model - Point estimate")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(prev_year = factor(year),
grp = "MAP Calibrated Model")) %>%
bind_rows(., dt_prev %>%
group_by(prev_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
prev_year = factor(prev_year))) %>%
bind_rows(., dt_prev %>%
group_by(prev_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
prev_year = factor(prev_year)))
ggplotly(ggplot(data = dt_prev,
aes(y = value, x = prev_year)) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = wprop_opioids_rx_tbl,
aes(y = value, x = prev_year, color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))Code
ggplotly(ggplot(data = dt_prev,
aes(x = value)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.3, position = 'identity') +
geom_vline(data = wprop_opioids_rx_tbl, aes(xintercept = value, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Prevalence of rx opioids use") + ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~prev_year, scales = "free"))Code
p_prev_sir <- ggplot(data = dt_prev,
aes(y = value, x = prev_year)) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = wprop_opioids_rx_tbl,
aes(y = value, x = prev_year, color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_prev_sir.png"))11.2 Total Deaths
Code
dt_deaths <- as.data.frame(params_sir_deaths_uncalib) %>%
pivot_longer(1:5, names_to = "death_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_deaths) %>%
pivot_longer(1:5, names_to = "death_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
ovrall_deaths_target_uncalib_tbl <- deaths_target_uncalib_tbl %>%
filter(target == "Total deaths") %>%
rename(value = target_val) %>%
mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
ifelse(group == "Target", "Target", group))),
death_year = factor(year)) %>%
select(-c(group, year)) %>%
bind_rows(., dt_deaths %>%
group_by(death_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
death_year = factor(death_year))) %>%
bind_rows(., dt_deaths %>%
group_by(death_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
death_year = factor(death_year)))
ggplotly(ggplot(data = dt_deaths,
aes(y = value, x = factor(death_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = ovrall_deaths_target_uncalib_tbl,
aes(y = value, x = factor(death_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))Code
ggplotly(ggplot(data = dt_deaths,
aes(x = value/1000)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.4, position = 'identity') +
# scale_fill_manual(values=c("#69b3a2", "#404080")) +
geom_vline(data = ovrall_deaths_target_uncalib_tbl,
aes(xintercept = value/1000, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Total deaths (in thousands)")+ ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
xlim(NA, 320)+
facet_wrap(~death_year, scales = "free"))Code
p_deaths_sir <- ggplot(data = dt_deaths,
aes(y = value, x = factor(death_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = ovrall_deaths_target_uncalib_tbl,
aes(y = value, x = factor(death_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_deaths_sir.png"))11.4 Total OAT counts
Code
dt_oat <- as.data.frame(params_sir_oat_uncalib) %>%
pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>%
mutate(group = "Uncalibrated Sample") %>%
bind_rows(., as.data.frame(opt_params_sir_oat) %>%
pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>%
mutate(group = "SIR Calibrated Sample"))
oat_target_uncalib_tbl <- oat_target_uncalib_tbl %>%
rename(value = target_val) %>%
mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
ifelse(group == "Target", "Target", group))),
oat_year = factor(year)) %>%
select(-c(group, year)) %>%
bind_rows(., dt_oat %>%
group_by(oat_year) %>%
summarize(value = mean(value)) %>%
mutate(grp = "SIR Calibrated Model - Mean",
oat_year = factor(oat_year))) %>%
bind_rows(., dt_oat %>%
group_by(oat_year) %>%
summarize(value = median(value)) %>%
mutate(grp = "SIR Calibrated Model - Median",
oat_year = factor(oat_year)))
ggplotly(ggplot(data = dt_oat,
aes(y = value, x = factor(oat_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = oat_target_uncalib_tbl,
aes(y = value, x = factor(oat_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group))Code
ggplotly(ggplot(data = dt_oat,
aes(x = value)) +
geom_histogram(aes(fill = group), color="#e9ecef",
alpha=0.4, position = 'identity') +
geom_vline(data = oat_target_uncalib_tbl,
aes(xintercept = value, color = grp))+
theme(axis.text.x = element_text(angle = 45)) +
xlab("Total OAT")+ ylab("") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~oat_year, scales = "free"))Code
p_oat_sir <- ggplot(data = dt_oat,
aes(y = value, x = factor(oat_year))) +
geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
geom_point(data = oat_target_uncalib_tbl,
aes(y = value, x = factor(oat_year), color = grp))+
xlab("Year") +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "") +
facet_wrap(~group)
ggsave(here("04_report/p_oat_sir.png"))12 Appendix 6: CBA results - point estimates
Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death
Code
trace_tbl <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_1$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl$mod <- factor(trace_tbl$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
mod_all_int_1$extra_sev_bi)
v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
mod_all_int_1$extra_mod_bi)
v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
mod_all_int_1$extra_od_deaths)
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)
for (j in 1:length(mod_names)){
trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_od[j]
trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
trace_tbl$cycle_num == 180 &
trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}
trace_tbl$state <- factor(trace_tbl$state,
levels = v_state_names,
labels = v_state_names)Code
trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "No Interventions") %>%
bind_rows(., as_tibble(mod_nalx_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Naloxone")) %>%
bind_rows(., as_tibble(mod_ss_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Safer Supply")) %>%
bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "Prescription Guidelines")) %>%
bind_rows(., as_tibble(mod_all_int_2$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count") %>%
mutate(mod = "All Interventions"))
trace_tbl_2$mod <- factor(trace_tbl_2$mod,
levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))
v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
mod_all_int_2$extra_sev_bi)
v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
mod_all_int_2$extra_mod_bi)
v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
mod_all_int_2$extra_od_deaths)
for (j in 1:length(mod_names)){
trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]
trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
trace_tbl_2$cycle_num == 180 &
trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}
trace_tbl_2$state <- factor(trace_tbl_2$state,
levels = v_state_names,
labels = v_state_names)
trace_tbl <- trace_tbl %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., trace_tbl_2 %>%
mutate(scenario = "Scenario 2"))12.1 Cumulative OD Deaths
12.1.0.0.1 Cumulative OD Deaths - Scenario 1
Code
p1 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p1)12.1.0.0.2 Cumulative OD Deaths - Scenario 2
Code
p2 <- ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p2)12.1.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
geom_point(size = 1) +
ylab("Cumulative OD Deaths") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario) +
labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))12.2 New OD Deaths per year
Code
inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <- rep(NA, 14)
for (i in 1:14){
yr <- 2015 + i
inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase,
inci_od_deaths_nalox,
inci_od_deaths_ss,
inci_od_deaths_pg,
inci_od_deaths_all),
scenario = rep("Scenario 1", 14*5)) %>%
bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
year = c(rep(2016:2029, 5)),
inci_od_deaths = c(inci_od_deaths_basecase_2,
inci_od_deaths_nalox_2,
inci_od_deaths_ss_2,
inci_od_deaths_pg_2,
inci_od_deaths_all_2),
scenario = rep("Scenario 2", 14*5)))
saveRDS(inc_od_death_tbl, file = here("01_data/inc_od_death_tbl_mod.RDS"))
target_od_deaths <- calib_target_tbl %>%
filter(group == "total_od_deaths") %>%
select(year, target) %>%
rename(inci_od_deaths = target) %>%
mutate(intervention = "Target")
saveRDS(target_od_deaths, file = here("01_data/inc_od_death_tbl_target.RDS"))
inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
levels = mod_names, labels = mod_names)12.2.0.0.1 Scenario 1
Code
p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p3)12.2.0.0.2 Scenario 2
Code
p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2")
ggplotly(p4)12.2.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(inc_od_death_tbl,
aes(x = year, y = inci_od_deaths, color = intervention)) +
geom_line() +
geom_point(size = 1) +
geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
ylab("New OD Deaths per year") + xlab("Year") +
labs(colour = "Intervention") +
scale_color_brewer(palette = "Dark2") +
facet_wrap(~scenario))12.3 Cumulative Deaths
12.3.0.0.1 Scenario 1
Code
p5 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 1") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p5)12.3.0.0.2 Scenario 2
Code
p6 <- ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
filter(scenario == "Scenario 2") %>%
mutate(year = rep(c(2015:2029), length(mod_names))),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention")
ggplotly(p6)12.3.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state == "BO_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = rep(c(2015:2029), length(mod_names)*2)),
aes(x = year, y = count, colour = mod)) +
geom_line() +
ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
labs(colour = "Intervention") +
facet_wrap(~scenario))12.4 Trace Plots
Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
ggplotly(ggplot(trace_tbl %>%
filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>%
filter(grepl("BPO", state)),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
"BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
"BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))Code
plotly::ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl %>%
filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
"BR_OAT_MAINT", "BR_OAT_SS",
"BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal) +
# facet_wrap(~mod) +
facet_grid(rows = vars(mod), cols = vars(scenario),
switch = "y"))12.5 CBA results
12.5.0.1 Scenario 1
Code
cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)
cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))
deaths_eff_tbl <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
total_death_oddeaths <- trace_tbl %>%
filter(scenario == "Scenario 1") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))
od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])
saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
od_deaths_diff, od_deaths_diff_per,
deaths_diff, deaths_diff_per), here("01_data/point_estiamte.RDS"))
saveRDS(total_death_oddeaths, here("01_data/total_deaths_oddeaths.RDS"))
costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")
effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs = "Discounted Net Present Costs \n in Millions (%)",
deaths = "Deaths (%)",
oddeaths = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 102 (0.01%) | 396 (0.01%) | -6420 (-4.05%) |
Safer Supply | 4233 (0.53%) | -13956 (-0.31%) | -3138 (-1.98%) |
Prescription Guidelines | -26103 (-3.28%) | 15252 (0.34%) | -1764 (-1.11%) |
All Interventions | -24876 (-3.13%) | 14115 (0.31%) | -8543 (-5.39%) |
12.5.0.2 Scenario 2
Code
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)
cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))
deaths_eff_tbl_2 <- trace_tbl %>%
filter(scenario == "Scenario 2") %>%
filter(cycle_num == 180) %>%
select(-cycle_num) %>%
pivot_wider(names_from = mod, values_from = count) %>%
filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>%
mutate(nalx_effect = Naloxone - `No Interventions`,
ss_effect = `Safer Supply` - `No Interventions`,
pg_effect = `Prescription Guidelines` - `No Interventions`,
all_effect = `All Interventions` - `No Interventions`,
nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>%
select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))
od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])
costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")
effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)
effects_tbl_2 |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Mean Change Compared to No Intervention")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
costs_2 = "Discounted Net Present Costs \n in Millions (%)",
deaths_2 = "Deaths (%)",
oddeaths_2 = "Overdose-related Deaths (%)")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Mean Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths (%) | Overdose-related Deaths (%) |
Naloxone | 85 (0.01%) | 303 (0.01%) | -4861 (-3.9%) |
Safer Supply | 4047 (0.51%) | -12743 (-0.28%) | -2435 (-1.95%) |
Prescription Guidelines | -26098 (-3.29%) | 15308 (0.34%) | -1485 (-1.19%) |
All Interventions | -24882 (-3.13%) | 14105 (0.31%) | -6607 (-5.3%) |
Code
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
deaths_diff_per, od_deaths_diff_per) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
deaths_diff_per_2, od_deaths_diff_per_2) %>%
pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>%
mutate(scenario = "Scenario 2")) %>%
mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))12.5.0.2.1 Scenario 1
Code
p7 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p7)12.5.0.2.2 Scenario 2
Code
p8 <- ggplot(data = tbl_effect_plot1 %>%
filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_minimal()
ggplotly(p8)12.5.0.2.3 Both
Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","OD-releated Deaths"))+
theme_bw() +
facet_wrap(~scenario))Code
tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
od_deaths_diff_per) %>%
mutate(scenario = "Scenario 1") %>%
bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
od_deaths_diff_per_2) %>%
rename(cost_diff_per = "cost_diff_per_2",
od_deaths_diff_per = "od_deaths_diff_per_2") %>%
mutate(scenario = "Scenario 2"))12.5.0.2.4 Scenario 1
Code
p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p9)12.5.0.2.5 Scenario 2
Code
p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_minimal()
ggplotly(p10)12.5.0.2.6 Both
Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
geom_point(aes(color = Intervention, shape = scenario)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
xlab("Mean Change in Costs Compared to No Intervention (%)") +
ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
theme_bw() +
facet_wrap(~scenario))12.6 Secondary Outcomes
Code
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)
prev_fun <- function(mod, i, v_states){
prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("prop_",
str_sub(v_yrs_prev, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12,
ncol = yrs_prev,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev, 3, 4)))))
for (j in 1:yrs_prev) {
yr <- (v_yrs_prev[1] - 1) + j
if (length(v_states) > 1)
{
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
} else {
prop_uncalib[, j] <- assign(
paste0("prop_", str_sub(yr, 3, 4)),
mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
v_states]) /
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
}
tot_pop_uncalib_tbl[, j] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
# final table for monthly prevalence over 15 years
prop_uncalib_tbl <- prop_uncalib %>%
pivot_longer(1:15, names_to = "grp", values_to = "value") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:15, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop)) %>%
mutate(intervention = mod_names[i])
# table with weighted yearly prevalence
prop_uncalib_wei_mean <- prop_uncalib_tbl %>%
group_by(year) %>%
summarise(value = weighted.mean(value, tot_pop_val)) %>%
ungroup() %>%
mutate(intervention = mod_names[i])
return (list(prop_uncalib_tbl = prop_uncalib_tbl,
prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}12.6.1 Severe Brain injury
Code
noint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
prev_sevbi_mon_tbl <- noint_prev_sevbi_mon_tbl %>%
bind_rows(., nalx_prev_sevbi_mon_tbl, ss_prev_sevbi_mon_tbl,
pg_prev_sevbi_mon_tbl, allint_prev_sevbi_mon_tbl)
prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>%
bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)
prev_sevbi_mon_tbl$year_mon <- factor(prev_sevbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_SEVERE_BI") %>%
filter(scenario == "Scenario 1") %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_sevbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of severe brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.2 Prevalence of moderate Brain injury
Code
noint_prev_modbi_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_MOD_BI")$prop_uncalib_tbl
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
prev_modbi_mon_tbl <- noint_prev_modbi_mon_tbl %>%
bind_rows(., nalx_prev_modbi_mon_tbl, ss_prev_modbi_mon_tbl,
pg_prev_modbi_mon_tbl, allint_prev_modbi_mon_tbl)
prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>%
bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)
prev_modbi_mon_tbl$year_mon <- factor(prev_modbi_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_MOD_BI") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_modbi_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of moderate brain injury") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.3 Prevalence of substance use treatment
Code
v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
"BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
"BR_OAT_SS", "BR_SS")
noint_prev_tx_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_tbl
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_tx)$prop_uncalib_tbl
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_tx)$prop_uncalib_tbl
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_tx)$prop_uncalib_tbl
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_tx)$prop_uncalib_tbl
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_tx)$prop_uncalib_wei_mean
prev_tx_mon_tbl <- noint_prev_tx_mon_tbl %>%
bind_rows(., nalx_prev_tx_mon_tbl, ss_prev_tx_mon_tbl,
pg_prev_tx_mon_tbl, allint_prev_tx_mon_tbl)
prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>%
bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)
prev_tx_mon_tbl$year_mon <- factor(prev_tx_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_prev_tx) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_tx_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of substance use treatment") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.4 Prevalence of prescription opioid use
Code
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target) %>%
rename(value = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")
noint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
prev_rx_opd_mon_tbl <- noint_prev_rx_opd_mon_tbl %>%
bind_rows(., nalx_prev_rx_opd_mon_tbl, ss_prev_rx_opd_mon_tbl,
pg_prev_rx_opd_mon_tbl, allint_prev_rx_opd_mon_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value, year_mon) %>%
mutate(intervention = "Target"))
prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>%
bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, value) %>%
mutate(intervention = "Target"))
prev_rx_opd_mon_tbl$year_mon <- factor(prev_rx_opd_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rx_opd_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of prescription opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.5 Prevalence of illicit use
Code
v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_tbl
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicit)$prop_uncalib_tbl
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicit)$prop_uncalib_tbl
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicit)$prop_uncalib_tbl
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicit)$prop_uncalib_tbl
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicit)$prop_uncalib_wei_mean
prev_illicituse_mon_tbl <- noint_prev_illicituse_mon_tbl %>%
bind_rows(., nalx_prev_illicituse_mon_tbl, ss_prev_illicituse_mon_tbl,
pg_prev_illicituse_mon_tbl, allint_prev_illicituse_mon_tbl)
prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>%
bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)
prev_illicituse_mon_tbl$year_mon <- factor(prev_illicituse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicit) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_illicituse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid use") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.6 Prevalence of misuse
Code
noint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BPO_MISUSE")$prop_uncalib_tbl
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
prev_rxmisuse_mon_tbl <- noint_prev_rxmisuse_mon_tbl %>%
bind_rows(., nalx_prev_rxmisuse_mon_tbl, ss_prev_rxmisuse_mon_tbl,
pg_prev_rxmisuse_mon_tbl, allint_prev_rxmisuse_mon_tbl)
prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>%
bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)
prev_rxmisuse_mon_tbl$year_mon <- factor(prev_rxmisuse_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BPO_MISUSE") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rxmisuse_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid misuse") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.7 Prevalence of overdose
Code
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_tbl
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_od)$prop_uncalib_tbl
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_od)$prop_uncalib_tbl
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_od)$prop_uncalib_tbl
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_od)$prop_uncalib_tbl
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_od)$prop_uncalib_wei_mean
prev_od_mon_tbl <- noint_prev_od_mon_tbl %>%
bind_rows(., nalx_prev_od_mon_tbl, ss_prev_od_mon_tbl,
pg_prev_od_mon_tbl, allint_prev_od_mon_tbl)
prev_od_yr_tbl <- noint_prev_od_yr_tbl %>%
bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)
prev_od_mon_tbl$year_mon <- factor(prev_od_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_od) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_od_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.7.1 prevalence of illicit overdose
Code
v_states_illicitod <- c("BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_tbl
noint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
nalx_prev_illicitod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = v_states_illicitod)$prop_uncalib_tbl
nalx_prev_illicitod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
ss_prev_illicitod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = v_states_illicitod)$prop_uncalib_tbl
ss_prev_illicitod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
pg_prev_illicitod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = v_states_illicitod)$prop_uncalib_tbl
pg_prev_illicitod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
allint_prev_illicitod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = v_states_illicitod)$prop_uncalib_tbl
allint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = v_states_illicitod)$prop_uncalib_wei_mean
prev_illicitod_mon_tbl <- noint_prev_illicitod_mon_tbl %>%
bind_rows(., nalx_prev_illicitod_mon_tbl, ss_prev_illicitod_mon_tbl,
pg_prev_illicitod_mon_tbl, allint_prev_illicitod_mon_tbl)
prev_illicitod_yr_tbl <- noint_prev_illicitod_yr_tbl %>%
bind_rows(., nalx_prev_illicitod_yr_tbl, ss_prev_illicitod_yr_tbl,
pg_prev_illicitod_yr_tbl, allint_prev_illicitod_yr_tbl)
prev_illicitod_mon_tbl$year_mon <- factor(prev_illicitod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% v_states_illicitod) %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_illicitod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))12.6.7.2 prevalence of rx overdose
Code
noint_prev_rxod_mon_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_tbl
noint_prev_rxod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
nalx_prev_rxod_mon_tbl <- prev_fun(mod = mod_nalx_1, i = 2,
v_states = "BO_OD_RX")$prop_uncalib_tbl
nalx_prev_rxod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
ss_prev_rxod_mon_tbl <- prev_fun(mod = mod_ss_1, i = 3,
v_states = "BO_OD_RX")$prop_uncalib_tbl
ss_prev_rxod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
pg_prev_rxod_mon_tbl <- prev_fun(mod = mod_pres_guid_1, i = 4,
v_states = "BO_OD_RX")$prop_uncalib_tbl
pg_prev_rxod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
allint_prev_rxod_mon_tbl <- prev_fun(mod = mod_all_int_1, i = 5,
v_states = "BO_OD_RX")$prop_uncalib_tbl
allint_prev_rxod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
v_states = "BO_OD_RX")$prop_uncalib_wei_mean
prev_rxod_mon_tbl <- noint_prev_rxod_mon_tbl %>%
bind_rows(., nalx_prev_rxod_mon_tbl, ss_prev_rxod_mon_tbl,
pg_prev_rxod_mon_tbl, allint_prev_rxod_mon_tbl)
prev_rxod_yr_tbl <- noint_prev_rxod_yr_tbl %>%
bind_rows(., nalx_prev_rxod_yr_tbl, ss_prev_rxod_yr_tbl,
pg_prev_rxod_yr_tbl, allint_prev_rxod_yr_tbl)
prev_rxod_mon_tbl$year_mon <- factor(prev_rxod_mon_tbl$year_mon,
levels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"
), 5),
labels = rep(
paste(
rep(v_yrs_prev,
each = 12),
month.abb[1:12],
sep = "_"), 5))
ggplotly(ggplot(trace_tbl %>%
filter(state %in% "BO_OD_RX") %>%
filter(scenario == "Scenario 1") %>%
group_by(cycle_num, mod) %>%
summarize(count = sum(count)) %>%
ungroup() %>%
left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
aes(x = cycle_num, y = count, colour = mod, key = year)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
scale_color_manual(values = colors_pal, name = "Intervention"))Code
ggplotly(ggplot() +
geom_point(data = prev_rxod_yr_tbl,
aes(x = factor(year), y = value,
color = intervention, fill = intervention)) +
xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid-related overdose") +
scale_fill_discrete(name = "Intervention")+
scale_color_discrete(name = "Intervention")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))13 Appendix 7: PSA results
Code
library(tidyverse)
library(plotly)
library(flextable)
theme_set(theme_minimal())
library(here)
m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))
mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)13.1 Distribution of change in primary outcomes compared to No intervention
13.1.2 Overall deaths
Code
ggplotly(ggplot(data = diff_psa_dt %>%
filter(grp == "death"), aes(x = diff, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Change in Deaths Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff %>%
filter(grp == "death"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())13.1.3 Costs
Code
ggplotly(ggplot(data = diff_psa_dt %>%
filter(grp == "cost"), aes(x = diff, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Change in Costs Compared to No Intervention (in millions)") + ylab("") +
geom_vline(data = point_est_diff %>%
filter(grp == "cost"), aes(xintercept = diff/1000000, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())13.2 Distribution of percent change in primary outcomes compared to No intervention
Code
diff_per_psa_dt <- outcomes_psa_dt %>%
mutate(cost_diff_per_nalox = cost_diff_per_nalox,
cost_diff_per_ss = cost_diff_per_ss,
cost_diff_per_pg = cost_diff_per_pg,
cost_diff_per_all = cost_diff_per_all) %>%
select(cost_diff_per_nalox, cost_diff_per_ss, cost_diff_per_pg, cost_diff_per_all,
death_diff_per_nalox, death_diff_per_ss, death_diff_per_pg, death_diff_per_all,
oddeath_diff_per_nalox, oddeath_diff_per_ss, oddeath_diff_per_pg, oddeath_diff_per_all) %>%
pivot_longer(1:12, names_to = "group", values_to = "diff_per") %>%
separate(col = group, into = c("grp", "type1", "type2", "Intervention"), sep = "_") %>%
select(-type1, -type2) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
point_est_diff_per <- point_est %>%
select(cost_diff_per, deaths_diff_per, od_deaths_diff_per) %>%
cbind.data.frame(Intervention = mod_names[-1]) %>%
pivot_longer(1:3,names_to = "grp", values_to = "diff") %>%
separate(grp, into = c("grp", "type"), sep = "_") %>% select(-type) %>%
mutate(grp = ifelse(grp == "od", "oddeath",
ifelse(grp == "deaths", "death", grp)),
type = "Point Estimate")13.2.2 Overall deaths
Code
ggplotly(ggplot(data = diff_per_psa_dt %>%
filter(grp == "death"), aes(x = diff_per, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("% Change in Deaths Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff_per %>%
filter(grp == "death"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())13.2.3 Costs
Code
ggplotly(ggplot(data = diff_per_psa_dt %>%
filter(grp == "cost"), aes(x = diff_per, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("% Change in Costs Compared to No Intervention") + ylab("") +
geom_vline(data = point_est_diff_per %>%
filter(grp == "cost"), aes(xintercept = diff, linetype = type))+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal())13.3 Tables of effects of interventions
Code
mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
round(ci_cost_diff_nalox[1]/1000000, 0),", ",
round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
round(ci_cost_diff_ss[1]/1000000, 0),", ",
round(ci_cost_diff_ss[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
round(ci_cost_diff_pg[1]/1000000, 0),", ",
round(ci_cost_diff_pg[2]/1000000, 0), "]"),
paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
round(ci_cost_diff_all[1]/1000000, 0),", ",
round(ci_cost_diff_all[2]/1000000, 0), "]"))
mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
round(ci_death_diff_nalox[1], 0),", ",
round(ci_death_diff_nalox[2], 0), "]"),
paste0(round(ci_death_diff_ss[3], 0), " [",
round(ci_death_diff_ss[1], 0),", ",
round(ci_death_diff_ss[2], 0), "]"),
paste0(round(ci_death_diff_pg[3], 0), " [",
round(ci_death_diff_pg[1], 0),", ",
round(ci_death_diff_pg[2], 0), "]"),
paste0(round(ci_death_diff_all[3], 0), " [",
round(ci_death_diff_all[1], 0),", ",
round(ci_death_diff_all[2], 0), "]"))
mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
round(ci_oddeath_diff_nalox[1], 0),", ",
round(ci_oddeath_diff_nalox[2], 0), "]"),
paste0(round(ci_oddeath_diff_ss[3], 0), " [",
round(ci_oddeath_diff_ss[1], 0),", ",
round(ci_oddeath_diff_ss[2], 0), "]"),
paste0(round(ci_oddeath_diff_pg[3], 0), " [",
round(ci_oddeath_diff_pg[1], 0),", ",
round(ci_oddeath_diff_pg[2], 0), "]"),
paste0(round(ci_oddeath_diff_all[3], 0), " [",
round(ci_oddeath_diff_all[1], 0),", ",
round(ci_oddeath_diff_all[2], 0), "]"))
effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
mean_death_diff, mean_oddeath_diff)
effects_tbl |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
mean_death_diff = "Deaths ",
mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present | Deaths | Opioid-related Overdose Deaths |
Naloxone | 209 [107, 362] | 625 [414, 926] | -10089 [-15009, -6697] |
Safer Supply | 4216 [3651, 4821] | -14180 [-17532, -11074] | -3221 [-4105, -2529] |
Prescription Guidelines | -25455 [-31480, -19866] | 14417 [8610, 20565] | -5458 [-11045, -1624] |
All Interventions | -21058 [-27059, -15403] | 960 [-5687, 7753] | -18544 [-28496, -11309] |
Code
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
round(ci_cost_diff_per_nalox[1], 2),", ",
round(ci_cost_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_cost_diff_per_ss[3], 2), " [",
round(ci_cost_diff_per_ss[1], 2),", ",
round(ci_cost_diff_per_ss[2], 2), "]%"),
paste0(round(ci_cost_diff_per_pg[3], 2), " [",
round(ci_cost_diff_per_pg[1], 2),", ",
round(ci_cost_diff_per_pg[2], 2), "]%"),
paste0(round(ci_cost_diff_per_all[3], 2), " [",
round(ci_cost_diff_per_all[1], 2),", ",
round(ci_cost_diff_per_all[2], 2), "]%"))
mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
round(ci_death_diff_per_nalox[1], 2),", ",
round(ci_death_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_death_diff_per_ss[3], 2), " [",
round(ci_death_diff_per_ss[1], 2),", ",
round(ci_death_diff_per_ss[2], 2), "]%"),
paste0(round(ci_death_diff_per_pg[3], 2), " [",
round(ci_death_diff_per_pg[1], 2),", ",
round(ci_death_diff_per_pg[2], 2), "]%"),
paste0(round(ci_death_diff_per_all[3], 2), " [",
round(ci_death_diff_per_all[1], 2),", ",
round(ci_death_diff_per_all[2], 2), "]%"))
mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
probs = c(0.025, 0.975, 0.5))
mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
round(ci_oddeath_diff_per_nalox[1], 2),", ",
round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
round(ci_oddeath_diff_per_ss[1], 2),", ",
round(ci_oddeath_diff_per_ss[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
round(ci_oddeath_diff_per_pg[1], 2),", ",
round(ci_oddeath_diff_per_pg[2], 2), "]%"),
paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
round(ci_oddeath_diff_per_all[1], 2),", ",
round(ci_oddeath_diff_per_all[2], 2), "]%"))
effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
mean_death_diff_per, mean_oddeath_diff_per)
effects_tbl_per |>
flextable() |>
add_header_row(colwidths = c(1, 3),
values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
set_header_labels(values = list(
`mod_names[-1]` = "Interventions",
mean_cost_diff_per = "Discounted Net Present Costs",
mean_death_diff_per = "Deaths ",
mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
theme_vanilla() |>
set_table_properties(layout = "autofit") |>
align_text_col(align = "center", header = TRUE, footer = TRUE)Percent Change Compared to No Intervention | |||
|---|---|---|---|
Interventions | Discounted Net Present Costs | Deaths | Opioid-related Overdose Deaths |
Naloxone | 0.03 [0.01, 0.05]% | 0.01 [0.01, 0.02]% | -3.46 [-4.13, -3.13]% |
Safer Supply | 0.54 [0.46, 0.61]% | -0.31 [-0.38, -0.24]% | -1.11 [-1.97, -0.66]% |
Prescription Guidelines | -3.23 [-3.94, -2.55]% | 0.32 [0.19, 0.45]% | -1.81 [-2.76, -0.93]% |
All Interventions | -2.67 [-3.39, -1.98]% | 0.02 [-0.13, 0.17]% | -6.43 [-7.16, -5.58]% |
13.4 Plots of effects of interventions
Code
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all,
ci_death_diff_per_nalox, ci_death_diff_per_ss,
ci_death_diff_per_pg, ci_death_diff_per_all,
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3)) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
# geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15, 16, 17, 18)) +
# scale_shape_manual(values = c(15, 16)) +
geom_hline(yintercept = 0) +
ylab("Percent Change Compared to No Intervention (%)") +
xlab("") +
scale_x_discrete(name ="",
labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
theme_minimal())Code
tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
ci_cost_diff_per_pg, ci_cost_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)) %>%
full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>%
pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
pivot_wider(names_from = "value", values_from = "per_diff") %>%
separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>%
select(-c(type1, type2, type3, group)), by = "Intervention") %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions", Intervention)))))
tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
geom_point(aes(color = Intervention), size = 1) +
geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention), height = 0) +
scale_color_brewer(palette = "Dark2") +
# scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_hline(yintercept = 0, linewidth = 0.25) +
geom_vline(xintercept = 0, linewidth = 0.25) +
xlab("Change in Costs Compared to No Intervention (%)") +
ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
theme_minimal())13.5 Posterior distributions for total Deaths and total OD-deaths
Code
death_psa_dt <- outcomes_psa_dt %>%
select(tot_death_no, tot_death_nalox, tot_death_ss, tot_death_pg, tot_death_all,
tot_oddeath_no, tot_oddeath_nalox, tot_oddeath_ss, tot_oddeath_pg, tot_oddeath_all) %>%
pivot_longer(1:10, names_to = "group", values_to = "total") %>%
separate(col = group, into = c("type", "grp", "Intervention"), sep = "_") %>%
select(-type) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions",
ifelse(Intervention == "no", "No Interventions", Intervention))))))
death_psa_dt$Intervention <- factor(death_psa_dt$Intervention,
levels = mod_names, labels = mod_names)
total_death_oddeaths <- total_death_oddeaths %>%
pivot_longer(3:7, names_to = "Intervention", values_to = "total") %>%
mutate(grp = ifelse(state == "BO_OD_DEATH", "oddeath", "death"),
type = "Point Estimate") %>%
select(-state)
ggplot(data = death_psa_dt %>%
filter(grp == "oddeath"), aes(x = total, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Total Opioid-related Overdose Deaths after 15 years") + ylab("") +
geom_vline(data = total_death_oddeaths %>%
filter(grp == "oddeath"), aes(xintercept = total, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal()
Code
ggplot(data = death_psa_dt %>%
filter(grp == "death"), aes(x = total/1000, fill = Intervention, color = Intervention)) +
geom_density(alpha = 0.9)+
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Total Deaths after 15 years (in thousands)") + ylab("") +
geom_vline(data = total_death_oddeaths %>%
filter(grp == "death"), aes(xintercept = total/1000, linetype = type), linewidth = 0.5)+
labs(linetype = "") +
facet_wrap(~Intervention, scales = "free_y") +
theme_minimal()
13.6 New OD Deaths per year
Code
outcomes_psa_dt_inci <- outcomes_psa_dt %>%
select(35:104)
inci_od_dt <- data.frame(nm = rep(NA, length(outcomes_psa_dt_inci)),
med = rep(NA, length(outcomes_psa_dt_inci)),
lb = rep(NA, length(outcomes_psa_dt_inci)),
ub = rep(NA, length(outcomes_psa_dt_inci)))
for (i in 1:length(outcomes_psa_dt_inci)){
inci_od_dt$nm[i] <- names(outcomes_psa_dt_inci)[i]
inci_od_dt[i, 2:4] <- as.numeric(quantile(outcomes_psa_dt_inci[i, ], c(0.5, 0.025, 0.975)))
}
inc_od_death_tbl_mod <- readRDS(here("01_data/inc_od_death_tbl_mod.RDS")) %>%
filter(scenario == "Scenario 1") %>%
select(-scenario) %>%
rename(Intervention = intervention,
Year = year,
value = inci_od_deaths) %>%
mutate(grp = "Point Estimate")
inci_od_dt_plot <- inci_od_dt %>%
separate(nm, into = c("type", "type1", "Intervention", "Year"), sep = "_") %>%
select(-c(type, type1)) %>%
rename(value = med) %>%
mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
ifelse(Intervention == "ss", "Safer Supply",
ifelse(Intervention == "pg", "Prescription Guidelines",
ifelse(Intervention == "all", "All Interventions",
ifelse(Intervention == "no", "No Interventions", Intervention))))),
Year = c(rep(2016:2029, 5)),
grp = "PSA results") %>%
bind_rows(., inc_od_death_tbl_mod) %>%
mutate(newgrp = paste0(Intervention, "_", grp))
inci_od_dt_plot$Intervention <- factor(inci_od_dt_plot$Intervention, levels = mod_names, labels = mod_names)
inc_od_death_tbl_target <- readRDS(here("01_data/inc_od_death_tbl_target.RDS")) %>%
mutate(grp = "Target")
ggplotly(ggplot(data = inci_od_dt_plot, aes(x = Year, y = value, color = grp, fill = grp)) +
geom_line() +
geom_point(data = inc_od_death_tbl_target, aes(x = year, y = inci_od_deaths, color = grp, fill = grp))+
scale_color_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3,2)]) +
scale_fill_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3, 2)]) +
# scale_color_brewer(palette = "Dark2") +
# scale_fill_brewer(palette = "Dark2") +
geom_ribbon(aes(ymin = lb, ymax = ub), alpha = 0.5) +
facet_wrap(~Intervention, scales = "free") +
theme(legend.title = element_blank()))14 References
Alarid-Escudero, Fernando, Richard F. MacLehose, Yadira Peralta, Karen M. Kuntz, and Eva A. Enns. 2018. “Nonidentifiability in Model Calibration and Implications for Medical Decision Making.” 2018. https://journals-sagepub-com.proxy3.library.mcgill.ca/doi/full/10.1177/0272989X18792283.
Busse, Jason W., Samantha Craigie, David N. Juurlink, D. Norman Buckley, Li Wang, Rachel J. Couban, Thomas Agoritsas, et al. 2017. “Guideline for Opioid Therapy and Chronic Noncancer Pain.” CMAJ 189 (18): E659–66. https://doi.org/10.1503/cmaj.170363.
Canada, Health. 2017. “Funding: Canadian Drugs and Substances Strategy.” Government of Canada. October 20, 2017. https://www.canada.ca/en/health-canada/services/substance-use/canadian-drugs-substances-strategy/funding.html.
———. 2020. “Federal Actions on Opioids to Date.” Transparency - other. Government of Canada. September 18, 2020. https://www.canada.ca/en/health-canada/services/opioids/federal-actions/overview.html.
———. 2021a. “Government of Canada Announces $20 Million to Help Communities Respond to Increasing Opioid-Related Overdoses.” News releases. Government of Canada. March 29, 2021. https://www.canada.ca/en/health-canada/news/2021/03/government-of-canada-announces-20-million-to-help-communities-respond-to-increasing-opioid-related-overdoses.html.
———. 2021b. “Safer Supply: Prescribed Medications as a Safer Alternative to Toxic Illegal Drugs.” Service description. July 22, 2021. https://www.canada.ca/en/health-canada/services/opioids/responding-canada-opioid-crisis/safer-supply.html.
Government of Canada, Statistics Canada. 2019. “Projected Population, by Projection Scenario, Age and Sex, as of July 1.” September 17, 2019. https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1710005701.
Health Canada, Health. 2020. “Canada’s Opioid Crisis (Fact Sheet).” Education and awareness. March 19, 2020. https://www.canada.ca/en/health-canada/services/publications/healthy-living/canada-opioid-crisis-fact-sheet.html.
Wyton, Moira. 2022. “Experts Reject BC’s Safe Supply Claims.” The Tyee; The Tyee. March 3, 2022. https://thetyee.ca/News/2022/03/03/Experts-Reject-BC-Safe-Supply-Claims/.